
library(mediation)
library(readstata13)
library(tidyverse)
library(ggplot2)
library(stargazer)
library(ggeffects)
#library(causalweight)
library(ivreg)
library(grid)
library(gridExtra)
library(interplot)
library(margins)
library(marginaleffects)
library(estimatr)
library(hettx)
library(sjPlot)

rm(list=ls(all=TRUE)) 
#setwd("/Users/hanilchang/Library/CloudStorage/Dropbox/My Working Paper/BLM and Covid19/Replication data")

setwd("D:/Dropbox/My Working Paper/BLM and Covid19/Replication data")

###########################
###########################
####        USA        ####
###########################
###########################
data.usa <- read.csv("replication_data_usa.csv")

table(data.usa$main.manipulation)
table(data.usa$main.manipulation, data.usa$democrat)
table(data.usa$main.manipulation, data.usa$republican)
t.test(manipulation.check ~ main.manipulation, data=data.usa)

## randomization check for partisans
data.usa <- subset(data.usa, republican == 1 | democrat ==1)

data.usa$diff.abs.ideo.my.party <- NA
data.usa$diff.abs.ideo.other.party <- NA


for (i in 1:nrow(data.usa)){
  if (data.usa$republican[i] == 1){
    data.usa$diff.abs.ideo.my.party[i] <- abs(data.usa$ideology.republican[i] - data.usa$conservative[i])
    data.usa$diff.abs.ideo.other.party[i] <- abs(data.usa$ideology.democratic[i] - data.usa$conservative[i])
  }
  if (data.usa$democrat [i] == 1){
    data.usa$diff.abs.ideo.my.party[i] <- abs(data.usa$ideology.democratic[i] - data.usa$conservative[i]) 
    data.usa$diff.abs.ideo.other.party[i] <- abs(data.usa$ideology.republican[i] - data.usa$conservative[i])
  }
}

table(data.usa$diff.abs.ideo.my.party)
table(data.usa$diff.abs.ideo.other.party)

t.test(age ~ main.manipulation, data=data.usa)
t.test(female ~ main.manipulation, data=data.usa)
t.test(income ~ main.manipulation, data=data.usa)
t.test(education ~ main.manipulation, data=data.usa)
t.test(white ~ main.manipulation, data=data.usa)
t.test(black ~ main.manipulation, data=data.usa)
t.test(asian ~ main.manipulation, data=data.usa)
t.test(hispanic ~ main.manipulation, data=data.usa)
t.test(covid.me.threatened ~ main.manipulation, data=data.usa)
t.test(republican ~ main.manipulation, data=data.usa)
t.test(diff.abs.ideo.my.party ~ main.manipulation, data=data.usa)
t.test(diff.abs.ideo.other.party ~ main.manipulation, data=data.usa)
t.test(poli.interest ~ main.manipulation, data=data.usa)
t.test(poli.participation ~ main.manipulation, data=data.usa)
t.test(poli.knowledge ~ main.manipulation, data=data.usa)
t.test(approval.trump ~ main.manipulation, data=data.usa)

##################
#### Table 1 ####
##################
med.fit.nocontrol.1 <- lm(diff.affect ~ main.manipulation,
              data = data.usa)
summary(med.fit.nocontrol.1)


out.fit.1 <- lm(social.bias ~ diff.affect + 
                  main.manipulation + republican + age + female + income + education + white + black + asian + hispanic + 
                  bea.region2 + bea.region3 + bea.region4 + bea.region5 + bea.region6 + bea.region7 + bea.region8 +
                  diff.abs.ideo.my.party + diff.abs.ideo.other.party +
                  poli.interest + 
                  poli.participation + 
                  poli.knowledge + 
                  approval.trump + 
                  covid.me.threatened,
                data = data.usa)
summary(out.fit.1)

predict.out.fit.1<- ggpredict(out.fit.1, terms=c("diff.affect"))
predict.out.fit.1

out.fit.2 <- lm(relative.stereotype ~ diff.affect + 
                  main.manipulation + republican + age + female + income + education + white + black + asian + hispanic + 
                  bea.region2 + bea.region3 + bea.region4 + bea.region5 + bea.region6 + bea.region7 + bea.region8 +
                  diff.abs.ideo.my.party + diff.abs.ideo.other.party +
                  poli.interest + 
                  poli.participation + 
                  poli.knowledge + 
                  approval.trump + 
                  covid.me.threatened,
                data = data.usa)
summary(out.fit.2)

predict.out.fit.2<- ggpredict(out.fit.2, terms=c("diff.affect"))
predict.out.fit.2


##################
### Appendix E ###
##################
out.fit.1.E <- lm(social.bias ~ 
                    main.manipulation + republican + age + female + income + education + white + black + asian + hispanic + 
                    bea.region2 + bea.region3 + bea.region4 + bea.region5 + bea.region6 + bea.region7 + bea.region8 +
                    diff.abs.ideo.my.party + diff.abs.ideo.other.party +
                    poli.interest + 
                    poli.participation + 
                    poli.knowledge + 
                    approval.trump + 
                    covid.me.threatened,
                  data = data.usa)
summary(out.fit.1.E)

out.fit.2.E <- lm(relative.stereotype ~ 
                    main.manipulation + republican + age + female + income + education + white + black + asian + hispanic + 
                    bea.region2 + bea.region3 + bea.region4 + bea.region5 + bea.region6 + bea.region7 + bea.region8 +
                    diff.abs.ideo.my.party + diff.abs.ideo.other.party +
                    poli.interest + 
                    poli.participation + 
                    poli.knowledge + 
                    approval.trump + 
                    covid.me.threatened,
                  data = data.usa)
summary(out.fit.2.E)



###############################
## Figure 1 and Appendix F ####
###############################
data.usa.sub <- dplyr::select(data.usa, 
                              social.bias,
                              relative.stereotype, 
                              diff.affect, 
                              main.manipulation, 
                              republican, 
                              age, female, income, education, white, black, asian, hispanic,
                              bea.region1, bea.region2,bea.region3,bea.region4,bea.region5,bea.region6,bea.region7,bea.region8,
                              diff.abs.ideo.my.party, 
                              diff.abs.ideo.other.party,
                              poli.interest, 
                              poli.participation, 
                              poli.knowledge, 
                              approval.trump,
                              manipulation.check, 
                              covid.me.threatened
) 
data.usa.sub <- na.omit(data.usa.sub)

set.seed(2020)

my.mediate <- function(x){
  med.fit <- lm(diff.affect ~  main.manipulation,
                data = data.usa.sub)
  
  out.fit <- lm(x ~ diff.affect +  main.manipulation + republican + age + female + income + education + white + black + asian + hispanic + 
                  bea.region2 + bea.region3 + bea.region4 + bea.region5 + bea.region6 + bea.region7 + bea.region8 +
                  diff.abs.ideo.my.party + diff.abs.ideo.other.party +
                  poli.interest + 
                  poli.participation + 
                  poli.knowledge + 
                  approval.trump + 
                  covid.me.threatened,
                data = data.usa.sub)
  
  med.out <- mediate(med.fit, out.fit, treat = "main.manipulation", mediator = "diff.affect", sim = 1000, robustSE = TRUE)
  out.list <- list(summary(med.fit), summary(out.fit), summary(med.out))
  return(out.list)
}

results <- lapply(data.usa.sub[, c("social.bias", "relative.stereotype")], function(x) my.mediate(x))
results ## Appendix F
 
comb.coef <- function(dv, model) {
  d <- data.frame(var = "Average Causal\nMediation Effect\n(ACME)",  coef = model$d1, l95 = model$d1.ci[1], u95 = model$d1.ci[2], pvalue = model$d1.p)
  z <- data.frame(var = "Average Direct\nEffect (ADE)",   coef = model$z1, l95 = model$z1.ci[1], u95 = model$z1.ci[2], pvalue = model$z1.p)
  t <- data.frame(var = "Total Effect", coef = model$tau.coef, l95 = model$tau.ci[1], u95 = model$tau.ci[2], pvalue = model$tau.p)
  n <- data.frame(var = "Explained", coef = model$n1, l95=model$n1.ci[1], u95 = model$n1.ci[2], pvalue = model$n1.p)
  out <- bind_rows(d, z, t, n)
  out$dv <- dv
  return(out)
}

out1 <- comb.coef("1. Social bias", results[["social.bias"]][3][[1]])  #ignore warnings 
out2 <- comb.coef("2. Relative stereotype", results[["relative.stereotype"]][3][[1]])  #ignore warnings 

comb.results <- bind_rows(out1, out2)
comb.results <- dplyr::filter(comb.results, var != "Explained")
comb.results

comb.plot.all <- ggplot(data = comb.results, aes(x=var, y = coef, fill = var)) + 
  geom_point(size = 2) + 
  geom_errorbar(aes(ymax = u95, ymin = l95), width = 0.1, size = 0.3) + 
  geom_hline(aes(yintercept = 0), size= 0.1, linetype =  2, color = "black") +
  theme_set(theme_bw()) + xlab("") + ylab("") + 
  facet_wrap( ~ dv, ncol=3, scales = "free_y") +
  theme(legend.position = "none") 
comb.plot.all 

a<-comb.plot.all 
ggsave("usa_mediation.png", a, width=7, height=3, scale=1.1)



###########################
###########################
###.   South Korea    ####
###########################
###########################
data.korea <- read.csv("replication_data_korea.csv")
data.korea <- subset(data.korea, partisanship.double.minjoo == 1 | partisanship.kookhim ==1)

data.korea$diff.abs.ideo.my.party <- NA
data.korea$diff.abs.ideo.other.party <- NA

for (i in 1:nrow(data.korea)){
  if (data.korea$partisanship.kookhim[i] == 1){
    data.korea$diff.abs.ideo.my.party[i] <- abs(data.korea$perceived.id.kookhim[i] - data.korea$self.id[i])
    data.korea$diff.abs.ideo.other.party[i] <- abs(data.korea$perceived.id.double.minjoo[i] - data.korea$self.id[i])
  }
  if (data.korea$partisanship.double.minjoo[i] == 1){
    data.korea$diff.abs.ideo.my.party[i] <- abs(data.korea$perceived.id.double.minjoo[i] - data.korea$self.id[i]) 
    data.korea$diff.abs.ideo.other.party[i] <- abs(data.korea$perceived.id.kookhim[i] - data.korea$self.id[i])
  }
}

table(data.korea$diff.abs.ideo.my.party)
table(data.korea$diff.abs.ideo.other.party)

t.test(age.cont~pol.exp.condition, data=data.korea)
t.test(female~pol.exp.condition, data=data.korea)
t.test(household.monthly.income~pol.exp.condition, data=data.korea)
t.test(education~pol.exp.condition, data=data.korea)
t.test(covid.anyone~pol.exp.condition, data=data.korea)
t.test(partisanship.kookhim~pol.exp.condition, data=data.korea)
t.test(diff.abs.ideo.my.party~pol.exp.condition, data=data.korea)
t.test(diff.abs.ideo.other.party~pol.exp.condition, data=data.korea)
t.test(poli.interest~pol.exp.condition, data=data.korea)
t.test(turnout.2017~pol.exp.condition, data=data.korea)
t.test(poli.knowledge.total~pol.exp.condition, data=data.korea)
t.test(gov.manage.covid ~pol.exp.condition, data=data.korea)


#################
### Table 2 ####
#################
med.fit.nocontrol.1 <- lm(diff.affect ~ pol.exp.condition,
                          data = data.korea)
summary(med.fit.nocontrol.1)

out.fit.1 <- lm(social.bias ~ diff.affect + pol.exp.condition + 
                  partisanship.kookhim + 
                  age.cont +  female +  household.monthly.income +  education +  
                  poli.interest +  
                  poli.knowledge.total + 
                  turnout.2017 + 
                  diff.abs.ideo.my.party + diff.abs.ideo.other.party + 
                  sudoguan +  kyeongsang +  jeolla +  choongcheong +  gangwon.jeju +  
                  gov.manage.covid + covid.anyone,
                data = data.korea)
summary(out.fit.1)

predict.out.fit.1<- ggpredict(out.fit.1, terms=c("diff.affect"))
predict.out.fit.1

out.fit.2 <- lm(relative.stereotype ~ diff.affect + pol.exp.condition + 
                  partisanship.kookhim + 
                  age.cont +  female +  household.monthly.income +  education + 
                  poli.interest +  
                  poli.knowledge.total + 
                  turnout.2017 + 
                  diff.abs.ideo.my.party + diff.abs.ideo.other.party + 
                  sudoguan +  kyeongsang +  jeolla +  choongcheong +  gangwon.jeju +  
                  gov.manage.covid + covid.anyone,
                data = data.korea)
summary(out.fit.2)

predict.out.fit.2<- ggpredict(out.fit.2, terms=c("diff.affect"))
predict.out.fit.2



###################
### Appendix K ####
###################
out.fit.1.E <- lm(social.bias ~ pol.exp.condition + 
                    partisanship.kookhim + 
                    age.cont +  female +  household.monthly.income +  education + 
                    poli.interest +  
                    poli.knowledge.total + 
                    turnout.2017 + 
                    diff.abs.ideo.my.party + diff.abs.ideo.other.party + 
                    sudoguan +  kyeongsang +  jeolla +  choongcheong +  gangwon.jeju +  
                    gov.manage.covid + covid.anyone,
                  data = data.korea)
summary(out.fit.1.E)

out.fit.2.E <- lm(relative.stereotype ~ pol.exp.condition + 
                    partisanship.kookhim + 
                    age.cont +  female +  household.monthly.income +  education + 
                    poli.interest +  
                    poli.knowledge.total + 
                    turnout.2017 + 
                    diff.abs.ideo.my.party + diff.abs.ideo.other.party + 
                    sudoguan +  kyeongsang +  jeolla +  choongcheong +  gangwon.jeju +  
                    gov.manage.covid + covid.anyone,
                  data = data.korea)
summary(out.fit.2.E)




####################################
####   Figure 2 and Appendix L  ###
####################################
data.korea.sub <- dplyr::select(data.korea, 
                                social.bias, 
                                relative.stereotype, 
                                diff.affect, 
                                pol.exp.condition, 
                                partisanship.kookhim,
                                age.cont, female, household.monthly.income, education, 
                                poli.interest, 
                                poli.knowledge.total,
                                turnout.2017, 
                                diff.abs.ideo.my.party, diff.abs.ideo.other.party, 
                                manipulation.check, 
                                seoul, sudoguan, kyeongsang, jeolla, choongcheong, gangwon.jeju, 
                                gov.manage.covid, covid.anyone)
data.korea.sub <- na.omit(data.korea.sub)

set.seed(2020)

my.mediate <- function(x){
  med.fit <- lm(diff.affect ~ pol.exp.condition,
                data = data.korea.sub)
  
  out.fit <- lm(x ~ diff.affect + pol.exp.condition + 
                  partisanship.kookhim + 
                  age.cont +  female +  household.monthly.income +  education +  
                  poli.interest +  
                  poli.knowledge.total + 
                  turnout.2017 + 
                  diff.abs.ideo.my.party + diff.abs.ideo.other.party + 
                  sudoguan +  kyeongsang +  jeolla +  choongcheong +  gangwon.jeju +  
                  gov.manage.covid +  covid.anyone,
                data = data.korea.sub)
  
  med.out <- mediate(med.fit, out.fit, treat = "pol.exp.condition", mediator = "diff.affect", sim = 1000, robustSE = TRUE)
  
  out.list <- list(summary(med.fit), summary(out.fit), summary(med.out))
  
  return(out.list)
}

results <- lapply(data.korea.sub[, c("social.bias", "relative.stereotype")], function(x) my.mediate(x))
results ## Appendix L

comb.coef <- function(dv, model) {
  d <- data.frame(var = "Average Causal\nMediation Effect\n(ACME)",  coef = model$d1, l95 = model$d1.ci[1], u95 = model$d1.ci[2], pvalue = model$d1.p)
  z <- data.frame(var = "Average Direct\nEffect (ADE)",   coef = model$z1, l95 = model$z1.ci[1], u95 = model$z1.ci[2], pvalue = model$z1.p)
  t <- data.frame(var = "Total Effect", coef = model$tau.coef, l95 = model$tau.ci[1], u95 = model$tau.ci[2], pvalue = model$tau.p)
  n <- data.frame(var = "Explained", coef = model$n1, l95=model$n1.ci[1], u95 = model$n1.ci[2], pvalue = model$n1.p)
  out <- bind_rows(d, z, t, n)
  out$dv <- dv
  return(out)
}

out1 <- comb.coef("1. Social bias", results[["social.bias"]][3][[1]])  #ignore warnings 
out2 <- comb.coef("2. Relative stereotype", results[["relative.stereotype"]][3][[1]])  #ignore warnings 

comb.results <- bind_rows(out1, out2)
comb.results <- dplyr::filter(comb.results, var != "Explained")

comb.plot.all <- ggplot(data = comb.results, aes(x=var, y = coef, fill = var)) + 
  geom_point(size = 2) + 
  geom_errorbar(aes(ymax = u95, ymin = l95), width = 0.1, size = 0.3) + 
  geom_hline(aes(yintercept = 0), size= 0.1, linetype =  2, color = "black") +
  theme_set(theme_bw()) + xlab("") + ylab("") + 
  facet_wrap( ~ dv, ncol=3, scales = "free_y") +
  theme(legend.position = "none")
comb.plot.all
a<-comb.plot.all 
ggsave("korea_mediation.png", a, width=7, height=3, scale=1.1)



#####################
###  Appendix M  ####
#####################
result.usa.control <- lm(diff.affect ~ main.manipulation + manipulation.check,
                         data = data.usa)
summary(result.usa.control)

result.korea.control <- lm(diff.affect ~ pol.exp.condition + manipulation.check,
                           data = data.korea)
summary(result.korea.control)



#####################
###  Appendix N  ####
#####################
usa.poli.interest <- lm(diff.affect ~ poli.interest*main.manipulation,
                        data = data.usa)
summary(usa.poli.interest)

usa.poli.knowledge <- lm(diff.affect ~ poli.knowledge*main.manipulation,
                         data = data.usa)
summary(usa.poli.knowledge)

korea.poli.interest <- lm(diff.affect ~ poli.interest*pol.exp.condition,
                          data = data.korea)
summary(korea.poli.interest)

korea.poli.knowledge.total <- lm(diff.affect ~ poli.knowledge.total*pol.exp.condition,
                                 data = data.korea)
summary(korea.poli.knowledge.total)



margin.plot.1<- plot_model(usa.poli.interest, type = "int", 
                           terms = c("main.manipulation[1,2]", "poli.interest [0, 1, 2, 3, 4]"), 
                           colors="bw",  ci.lvl= 0.95, title = "",
                           legend.title = "Experimental vignette", 
                           axis.title = "Predicted difference", 
                           show.legend = TRUE) + 
  coord_cartesian(expand = TRUE) + scale_x_continuous(limits = c(0, 4), breaks = c(0 ,1, 2, 3, 4)) + 
  xlab("Political interest") + 
  theme(legend.position=c(0.8,0.2)) + scale_linetype_discrete(name  ="Experimental vignette", 
                                                              breaks=c("1", "2"),
                                                              labels=c("Moderate", "Extreme"))
margin.plot.1
ggsave("usa_poli_interest.png", margin.plot.1, width=6, height=3, scale=1.1)

margin.plot.2<- plot_model(usa.poli.knowledge, type = "int", 
                           terms = c("main.manipulation[1,2]", "poli.knowledge [0, 1, 2, 3, 4, 5]"), 
                           colors="bw",  ci.lvl= 0.95, title = "",
                           legend.title = "Experimental vignette", 
                           axis.title = "Predicted difference", 
                           show.legend = TRUE) + 
  coord_cartesian(expand = TRUE)  + scale_x_continuous(limits = c(0, 5), breaks = c(0 ,1, 2, 3, 4, 5)) + 
  xlab("Political knowldege") + 
  theme(legend.position=c(0.8,0.2)) + scale_linetype_discrete(name  ="Experimental vignette", 
                                                              breaks=c("1", "2"),
                                                              labels=c("Moderate", "Extreme"))
margin.plot.2
ggsave("usa_poli_knowledge.png", margin.plot.2, width=6, height=3, scale=1.1)


margin.plot.1<- plot_model(korea.poli.interest, type = "int", 
                           terms = c("pol.exp.condition[1,2]", "poli.interest [0, 1, 2, 3]"), 
                           colors="bw",  ci.lvl= 0.95, title = "",
                           legend.title = "Experimental vignette", 
                           axis.title = "Predicted difference", 
                           show.legend = TRUE) + 
  coord_cartesian(expand = TRUE) + scale_x_continuous(limits = c(0, 3), breaks = c(0 ,1, 2, 3)) + 
  xlab("Political interest") + 
  theme(legend.position=c(0.8,0.2)) + scale_linetype_discrete(name  ="Experimental vignette", 
                                                              breaks=c("1", "2"),
                                                              labels=c("Moderate", "Extreme"))
margin.plot.1
ggsave("korea_poli_interest.png", margin.plot.1, width=6, height=3, scale=1.1)

margin.plot.2<- plot_model(korea.poli.knowledge.total, type = "int", 
                           terms = c("pol.exp.condition[1,2]", "poli.knowledge.total [0, 1, 2, 3, 4, 5, 6]"), 
                           colors="bw",  ci.lvl= 0.95, title = "",
                           legend.title = "Experimental vignette", 
                           axis.title = "Predicted difference", 
                           show.legend = TRUE) + 
  coord_cartesian(expand = TRUE) + scale_x_continuous(limits = c(0, 6), breaks = c(0 ,1, 2, 3, 4, 5, 6)) + 
  xlab("Political knowledge") + 
  theme(legend.position=c(0.8,0.2)) + scale_linetype_discrete(name  ="Experimental vignette", 
                                                              breaks=c("1", "2"),
                                                              labels=c("Moderate", "Extreme"))
margin.plot.2
ggsave("korea_poli_knowledge.png", margin.plot.2, width=6, height=3, scale=1.1)


#####################
###  Appendix O  ####
#####################
data.usa.sub <- dplyr::select(data.usa, 
                              social.bias,
                              relative.stereotype, 
                              diff.affect, 
                              main.manipulation, 
                              republican, 
                              age, female, income, education, white, black, asian, hispanic,
                              bea.region1, bea.region2,bea.region3,bea.region4,bea.region5,bea.region6,bea.region7,bea.region8,
                              diff.abs.ideo.my.party, 
                              diff.abs.ideo.other.party,
                              poli.interest, 
                              poli.participation, 
                              poli.knowledge, 
                              approval.trump,
                              covid.me.threatened) 

data.korea.sub <- dplyr::select(data.korea, 
                                social.bias, 
                                relative.stereotype, 
                                diff.affect, 
                                pol.exp.condition, 
                                partisanship.kookhim,
                                age.cont, female, household.monthly.income, education, 
                                poli.interest, 
                                poli.knowledge.total,
                                turnout.2017, 
                                diff.abs.ideo.my.party, diff.abs.ideo.other.party, 
                                seoul, sudoguan, kyeongsang, jeolla, choongcheong, gangwon.jeju, 
                                gov.manage.covid, covid.anyone)

### normalize
scale.data.usa <- as.data.frame(scale(data.usa.sub))
scale.data.usa$country <- 1
scale.data.usa$main.manipulation <- data.usa.sub$main.manipulation
summary(scale.data.usa)

scale.data.korea <- as.data.frame(scale(data.korea.sub))
scale.data.korea$country <- 0
scale.data.korea$pol.exp.condition <- data.korea.sub$pol.exp.condition
summary(scale.data.korea)

scale.data.usa.sub <- dplyr::select(scale.data.usa, 
                                    social.bias,
                                    relative.stereotype, 
                                    diff.affect, 
                                    main.manipulation, 
                                    republican, 
                                    age, female, income, education, 
                                    diff.abs.ideo.my.party, 
                                    diff.abs.ideo.other.party,
                                    poli.interest, 
                                    poli.participation, 
                                    poli.knowledge, 
                                    approval.trump,
                                    country)

scale.data.korea.sub <- dplyr::select(scale.data.korea, 
                                      social.bias, 
                                      relative.stereotype, 
                                      diff.affect, 
                                      pol.exp.condition, 
                                      partisanship.kookhim,
                                      age.cont, female, household.monthly.income, education, 
                                      diff.abs.ideo.my.party, 
                                      diff.abs.ideo.other.party, 
                                      poli.interest, 
                                      turnout.2017, 
                                      poli.knowledge.total,
                                      gov.manage.covid, 
                                      country)

scale.data.usa.sub <- rename(scale.data.usa.sub, 
                             incumbency.approval = "approval.trump")
colnames(scale.data.usa.sub)

scale.data.korea.sub <- rename(scale.data.korea.sub, 
                               social.bias = "social.bias", 
                               main.manipulation = "pol.exp.condition", 
                               republican = "partisanship.kookhim", 
                               age = "age.cont", 
                               income = "household.monthly.income", 
                               poli.interest = "poli.interest", 
                               poli.participation = "turnout.2017",
                               poli.knowledge = "poli.knowledge.total", 
                               incumbency.approval = "gov.manage.covid")

colnames(scale.data.korea.sub)
data.both <- rbind(scale.data.usa.sub, scale.data.korea.sub )
summary(data.both)
table(data.both$main.manipulation)
data.both.sub <- na.omit(data.both)


#### Table M-1
med.fit.nocontrol.1 <- lm(diff.affect ~ main.manipulation + country,
                          data = data.both)
summary(med.fit.nocontrol.1)

out.fit.1 <- lm(social.bias ~ diff.affect + 
                  main.manipulation + republican + age + female + income + education + 
                  diff.abs.ideo.my.party + diff.abs.ideo.other.party +
                  poli.interest + 
                  poli.participation + 
                  poli.knowledge + 
                  incumbency.approval + 
                  country,
                data = data.both)
summary(out.fit.1)

predict.out.fit.1<- ggpredict(out.fit.1, terms=c("diff.affect"))
predict.out.fit.1

out.fit.2 <- lm(relative.stereotype ~ diff.affect + 
                  main.manipulation + republican + age + female + income + education + 
                  diff.abs.ideo.my.party + diff.abs.ideo.other.party +
                  poli.interest + 
                  poli.participation + 
                  poli.knowledge + 
                  incumbency.approval + 
                  country,
                data = data.both)
summary(out.fit.2)

predict.out.fit.2<- ggpredict(out.fit.2, terms=c("diff.affect"))
predict.out.fit.2


### Table M-2 and Figure M-1
set.seed(2020)
my.mediate <- function(x){
  med.fit <- lm(diff.affect ~  main.manipulation + country,
                data = data.both.sub)
  
  out.fit <- lm(x ~ diff.affect +  main.manipulation + republican + age + female + income + education + 
                  diff.abs.ideo.my.party + diff.abs.ideo.other.party +
                  poli.interest + 
                  poli.participation + 
                  poli.knowledge + 
                  incumbency.approval + 
                  country
                ,
                data = data.both.sub)
  
  med.out <- mediate(med.fit, out.fit, treat = "main.manipulation", mediator = "diff.affect", sim = 1000, robustSE = TRUE)
  out.list <- list(summary(med.fit), summary(out.fit), summary(med.out))
  return(out.list)
}

results <- lapply(data.both.sub[, c("social.bias", "relative.stereotype")], function(x) my.mediate(x))

comb.coef <- function(dv, model) {
  d <- data.frame(var = "Average Causal\nMediation Effect\n(ACME)",  coef = model$d1, l95 = model$d1.ci[1], u95 = model$d1.ci[2], pvalue = model$d1.p)
  z <- data.frame(var = "Average Direct\nEffect (ADE)",   coef = model$z1, l95 = model$z1.ci[1], u95 = model$z1.ci[2], pvalue = model$z1.p)
  t <- data.frame(var = "Total Effect", coef = model$tau.coef, l95 = model$tau.ci[1], u95 = model$tau.ci[2], pvalue = model$tau.p)
  n <- data.frame(var = "Explained", coef = model$n1, l95=model$n1.ci[1], u95 = model$n1.ci[2], pvalue = model$n1.p)
  out <- bind_rows(d, z, t, n)
  out$dv <- dv
  return(out)
}

out1 <- comb.coef("1. Social bias", results[["social.bias"]][3][[1]])  #ignore warnings 
out2 <- comb.coef("2. Relative stereotype", results[["relative.stereotype"]][3][[1]])  #ignore warnings 

comb.results <- bind_rows(out1, out2)
comb.results <- dplyr::filter(comb.results, var != "Explained")
comb.results

comb.plot.all <- ggplot(data = comb.results, aes(x=var, y = coef, fill = var)) + 
  geom_point(size = 2) + 
  geom_errorbar(aes(ymax = u95, ymin = l95), width = 0.1, size = 0.3) + 
  geom_hline(aes(yintercept = 0), size= 0.1, linetype =  2, color = "black") +
  theme_set(theme_bw()) + xlab("") + ylab("") + 
  facet_wrap( ~ dv, ncol=3, scales = "free_y") +
  theme(legend.position = "none") 
comb.plot.all 

a<- comb.plot.all 
ggsave("meta_analysis_mediation.png", a, width=7, height=3, scale=1.1)



